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INTRODUCTION 

Warheads  impacting  ships,  concrete  bunkers  or  other  targets  are  often  required 
to  survive  perforation  or  penetration  in  detonable  condition  to  successfully  complete 
their  missions.  The  design  of  such  warheads  requires  a  knowledge  of  the  stresses  and 
deformations  induced  by  impact  from  which  an  assessment  of  survivability  can  be 
made.  Finite  element  structural  analysis  procedures  offer  an  attractive  way  of  deter¬ 
mining  these  quantities. 

Analyses  of  warheads  impacting  targets  at  normal  obliquity*  have  been  per¬ 
formed  at  the  Naval  Weapons  Center  using  the  two-dimensional  finite  element  code, 
HONDO  II.2  Analyses  of  warheads  impacting  targets  at  non-normal  obliquity  require  a 
full  three-dimensional  treatment  and  can  be  extremely  expensive  using  currently  avail¬ 
able  codes. 

This  report  documents  the  development  of  a  three-dimensional  finite  element 
code  for  the  analysis  of  axisymmetric  warheads  impacting  targets  at  obliquity.  The  aim 
of  this  effort  was  to  produce  a  specialized  and  efficient  code  well-adapted  to  the  Naval 
Weapons  Center  computer  system  (UNIVAC  1100/82)  and  capable  of  providing  useful 
information  to  the  warhead  designer  at  a  reasonable  cost.  In  addition,  it  was  hoped  to 
make  the  new  code  as  compatible  as  possible  with  locally  written  pre-  and  post¬ 
processors  for  the  HONDO  II  program  to  minimize  the  extra  effort  required  for 
oblique  impact  analyses. 

For  presentation  purposes  this  work  can  be  divided  into  three  parts: 

1.  Finite  element  theory. 

2.  Computer  implementation  of  theory;  i.e.,  writing  and  checkout  of  code. 

3.  Employment  of  the  code  in  the  solution  of  a  realistic  projectile  impact 
problem. 


*  J.  C.  Schulz.  “Finite  Element  Analysis  of  a  Kinetic  Energy  Warhead  Penetrating  Concrete.”  in  Proceed¬ 
ings  of  the  4th  International  Sym/x>sium  on  Ballistics,  Monterey,  California,  October  1978.  Washington,  D.C.,  Amer¬ 
ican  Defense  Preparedness  Association,  1978.  (Publication  UNCLASSIFIED.) 

2  Sandia  Laboratories.  HONDO  II.  A  Finite  Element  Computer  Program  for  large  Deformation  Dynamic 
Response  of  Axisymmetric  Solids,  by  S.  W.  Key,  Z.  E.  Bcisinger  and  R.  D.  Krieg.  Albuquerque,  N.  Mex.,  October 
1978.  (SAND78-0422,  publication  UNCLASSIFIED.) 
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These  are  discussed  in  subsequent  sections.  User  instructions  and  other  pertinent 
information  are  given  in  appendixes. 
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FINITE  ELEMENT  THEORY 


PROBLEM  DESCRIPTION 

A  solid  body  of  revolution  generated  by  rotating  a  planar  area  around  an  axis 
is  shown  in  Figure  1 .  Points  in  the  body  may  be  specified  by  rectangular  cartesian 
coordinates  (x^  i  =  1,  2,  3)  or  cylindrical  coordinates  (r,  6.  7).  Cylindrical  coordinates 
will  be  used  for  convenience  in  data  input,  while  rectangular  cartesian  coordinates  will 
be  used  for  the  finite  element  solution. 

At  time  t  =  0.  the  initial  coordinates  of  points  in  the  body  are  (Xs)  and  the 
body  is  assumed  to  be  strain-free.  Subsequently,  a  normal  pressure  loading  (not  neces¬ 
sarily  axisymmetric)  may  be  applied  resulting  in  traction  components  (S,)  acting  over 
a  portion  of  the  surface,  As.  The  body  may  also  have  an  initial  homogeneous  motion 
such  that  all  points  have  axial  velocity  components  (0,  0,  V3)  at  t  =  0.  Displacement 
boundary  conditions  (limited  to  the  setting  of  one  or  more  displacement  components 
equal  to  zero)  may  be  specified  over  some  reeion,  Ad. 


*2 


FIGURE  1.  Solid  Body  of  Revolution  Showing  Rectangular  Cartesian  and 
Cylindrical  Coordinate  Systems  and  Generating  Plane. 


4 


NWC  TP  6369 


As  a  result  of  the  pressure  loading  and  initial  velocity  and  subject  to  the  dis¬ 
placement  boundary  conditions,  the  body  will  undergo  a  motion 

xi  =  xi(X1,X2,X3,t)  (1) 

after  which  it  may  no  longer  be  axisymmetric.  The  problem  is  to  determine  the 
motion  and  to  ascertain  the  displacement  and  stresses  throughout  the  body  as  func¬ 
tions  of  time.  An  analytic  solution  is,  of  course,  not  feasible  except  in  certain  special 
cases,  and  some  form  of  approximate  solution  is  required.  The  finite  element  method 
is  used  to  obtain  such  a  solution. 


STARTING  POINT 

A  general  three-dimensional  theory  for  the  large  displacement,  inelastic  analysis 
of  solid  bodies  subjected  to  impulsive  loads  is  given  by  Key  and  then  specialized  to 
solids  of  revolution  subjected  to  axisymmetric  loads  to  form  the  basis  for  the  HONDO 
II  finite  element  code.2  The  same  general  theory  is  derived  here;  the  only  specializa¬ 
tion  being  to  a  rectangular  cartesian  coordinate  system. 


FINITE  ELEMENT  IDEALIZATION 

|  The  generating  planar  area  is  subdivided  into  quadrilateral  regions  and  then 

rotated  around  the  axis  to  form  a  set  of  annular  rings  of  quadrilateral  cross  section. 
The  rings  are  cut  by  a  series  of  radial  planes.  (Equal  spacing  of  the  planes  in  the  cir¬ 
cumferential  direction  is  used  in  the  computer  implementation,  although  not  required 
by  the  theory.)  The  body  is  thus  divided  into  hexahedral  elements  as  indicated  in 
I  Figure  2. 

A  typical  element,  as  shown  in  Figure  3,  can  be  specified  by  its  eight  comer 
nodal  points.  By  convention  the  numbering  must  be  counterclockwise  around  the 
quadrilateral  faces  in  consecutive  radial  planes.  A  local  coordinate  system  is  introduced 
so  that  the  initial  location  of  a  typical  point  inside  the  element  can  be  given  in  terms 
I  of  its  global  coordinates  (X|)  or  its  local  coordinates  (aj).  The  transformation  equations 

are 

8 

=  £  Xij0j(ai  ,a2,a3)  (2) 

I  j_1 

where  (xy)  are  the  initial  global  coordinates  of  the  nodes  (in  the  element  numbering 
system)  and  0j(ai,a2,a3)  are  shape  functions  defined  by 
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FIGURE  2.  Solid  of  Revolution  Partially  Cut  Away 
Showing  Division  Into  Finite  Elements. 


FIGURE  3.  Typical  Finite  Element  Showing  Element 
Node  Numbering  Convention. 
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0]  =  (1  -  a j )( 1  -  a2 )( 1  -  ^3 )/8 

02  ~  (1  +  a,  )(1  -  a2)(l  -  a3)/8 

03  ~  (1  +  a i )( 1  +  a2)(l  -  a3 )/8 
04  =  (1  -  a, XI  +  a2)(l  -  a3)/8 

05  =  (1  -  a, XI  ~  a2)(l  +  a3 )/8 

06  =  (1  +  a, XI  ~  a2)(l  +  a3)/8 

07  =  (1  +  a ! )( 1  +  a2)(I  +  a3)/8 

08  =  (1  -  a,  XI  +  a2Xl  +  a3)/8  (3) 

Subdivision  of  the  body  into  elements  and  the  introduction  of  local  coordinates 
are  matters  of  definition  and  involve  no  approximation  (except  at  the  surface  of  the 
body  where  the  planar  faces  of  the  elements  will  not  exactly  fit  the  curved  surface). 
At  this  point,  however,  the  fundamental  finite  element  idealization  is  made.  The 
displaced  positions  of  points  inside  the  deformed  element  are  approximated  by 


Xj  =  5Z  Xij0j  (4) 

j=l 

where  (x^)  are  the  global  coordinates  of  the  nodes.  With  this  approximation  the 
number  of  degrees  of  freedom  in  the  idealized  body  is  reduced  from  an  uncountable 
infinity  to  three  times  the  number  of  nodes  (unless  displacement  boundary  conditions 
are  imposed  at  some  of  the  nodes).  Displacement  compatibility  at  the  common  faces 
of  adjoining  elements  is  maintained. 

Equations  4  give  the  positions  of  points  within  any  one  element  in  terms  of 
the  positions  of  its  nodes.  A  different  set  of  equations  must  be  written  for  each 
element.  Equivalently,  a  single  set  of  equations  can  be  written  which  gives  the  posi¬ 
tions  of  points  throughout  the  entire  body.  Thus, 

Np 

x,  =  Xjj4>j  (5) 

j=l 

where  Np  is  the  total  number  of  nodal  points  in  the  body,  (Xjj)  are  the  coordinates 
of  the  jth  node  (in  a  global  node  numbering  system  rather  than  the  element  node 
numbering  system  of  Equations  4)  and  4>j  are  the  basis  functions  defined  by: 

4>j  =  0  in  elements  for  which  the  jth  node  is  not  a  comer  node, 

4>j  =  0k  elements  for  which  the  jth  node  is  a  comer  node,  where  k  is  the 
number  of  the  jth  node  in  the  element  node  numbering  system. 
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EQUATIONS  OF  MOTION 

Equations  of  motion  for  the  finite  element  model  of  the  body  can  be  obtained 
using  the  principle  of  virtual  displacements.  This  principle  requires  that  the  differential 
form 


r  ..  r  3(6xi)  r 

5tt=  /  px15xidV+  /  dV-  /  Sjfixj. 

J\j  J\J  OAm  •'A 


must  vanish  at  all  points  along  the  path  of  the  motion  for  all  virtual  displacements 
(5Xi>  satisfying  the  displacement  boundary  conditions.  The  virtual  displacements  can  be 
expressed  in  terms  of  the  virtual  displacements  of  the  nodes  as 

Np 

8xi  =  2iSxij*j 

j=l 

If  the  virtual  displacement  components  of  the  kth  node  are  in  turn  set  equal  to 
unity  while  all  other  nodal  virtual  displacement  components  are  set  equal  to  zero,  the 
following  equations  of  motion  result 


2$  r  r  r 

£  m  Jr^ a  "X  v  ■ +  X  s‘*kdA 


At  this  point  two  approximations  are  introduced  to  make  the  equations  of 
motion  more  manageable. 

First,  constant  stress  is  assumed  within  each  element.  This  approximation 
greatly  reduces  the  amount  of  storage  and  computational  time  required  for  solution; 
however,  it  can  lead  to  the  presence  of  “keystoning”  or  “zero-energy”  deformation 
modes.  These  can  be  suppressed  by  the  introduction  of  a  suitable  viscosity. 

Second,  an  approximate  diagonal  mass  matrix  is  formed  by  lumping  the  off- 
diagonal  terms  on  the  left-hand-side  of  Equations  8,  that  is,  the  mass  at  the  kth  node 
is  taken  as 


p 

=  Y,  f  p<M>kdV=  /  p4>kdV 

j=i  Jy  Jv 


where  the  second  equality  can  be  proved  by  integrating  over  individual  elements, 
expressing  the  basis  functions  in  terms  of  the  shape  functions  and  using  the  identity 
8 

E  +j  =  !• 

J=l 
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With  these  approximations,  Equations  8  become 


r  f.  34>k  i  /• 

myXjfc  =  Y1  |gim  J  ^r~dVj  +  J  Si<t>kdA  (k  not  summed)  (10) 

where  Ne  is  the  number  of  elements  and  the  subscript  j  refers  to  the  jth  element. 


TIME  INTEGRATION  SCHEME 

Equations  10  are  integrated  approximately  using  an  explicit  central  difference 
scheme.  It  is  assumed  that  the  nodal  coordinates  (x"k)  at  time  tn  are  known,  and  that 
the  nodal  velocity  components  (xfk”l,/2)  during  the  previous  time  step  Atn~'/~  are 
also  known.  The  nodal  acceleration  components  are  calculated  using  Equations  10. 
These  components  are  assumed  constant  over  the  new  time  step  Atn+1^2,  and  the 
nodal  velocity  components  are  updated  using 

(Atn_  1  /2  +  Atn+ 1  /21 

in-  1/2  =  1/2  +  x-  ^ >  (ID 

Similarly,  the  nodal  coordinates  are  updated  by 

xfk+l  =x['k  +  xink+I/2Atn+l/2  (12) 

Once  the  velocities  and  positions  of  nodal  points  have  been  updated,  the  strain 
rate  (more  properly  the  deformation  rate)  at  the  center  of  each  element  can  be  deter¬ 
mined  and  used  to  update  the  element  stresses  through  the  constitutive  relations.  Then 
the  procedure  can  be  repeated  for  the  next  time  step. 

To  damp  out  localized  oscillations  of  the  concentrated  masses  associated  with 
the  passage  of  shock  waves  through  the  material,  an  artificial  viscosity  pressure  is  intro¬ 
duced  in  the  form 


CiPWXi  +  c2P^ind2.ifdv<o 
0,  ifdv>  0 


Cj ,  Ct  =  constants 

jn  =  minimum  of  element  edge  and  diagonal  lengths 
c  =  uniaxial  sound  speed  in  element 

dv  =  d jj,  the  volumetric  deformation  rate  where  the  deformation  rate  compo¬ 
nents  (d(j)  will  be  defined  in  Equations  18. 
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The  quadratic  term  is  intended  to  damp  out  gradients  associated  with  shock  waves,  the 
linear  term  is  added  to  eliminate  additional  spurious  oscillations.  This  artificial  viscosity 
pressure  is  added  to  the  normal  stress  components  (»|(.  02  2-  °3s)  calculated  from  the 
material  constitutive  relations  before  introduction  into  Equations  10. 

This  integration  scheme  is  only  conditionally  stable.  To  ensure  stability  the 
time  step  must  be  kept  below  a  critical  value  dependent  on  the  material  sound  speed 
and  the  element  dimensions.  A  critical  time  step  for  each  element  is  computed  at  each 
time  step  using  the  formula 


C  |  C  +  C's^m m  klv  i  +  \/(Cj  C  +  Cs ^min  ldv  ^  +  V2 

which  accounts  for  artificial  viscosity  effects.  The  time  step  used  for  the  next  integra¬ 
tion  step  is  taken  as  0.9  times  the  minimum  of  the  critical  time  step  values  for  all  the 
elements. 


EVALUATION  OF  INTEGRALS 

The  volume  integrals  in  Equations  10  are  evaluated  approximately  within  each 
element  by  8-point  Gaussian  quadrature.  Thus 


(15) 


where 


X 
U I 


=  an  arbitrary  integrand 


i)Xj  3xj  0X|j 

eijk  3a  |  9aT  Saj' 


the  Jacobian  determinant  of  the  transformation  of 


Equations  4 

k  =  subscript  indicating  that  the  subscripted  quantities  are  to  be  evaluated 
at  the  Gaussian  quadrature  points,  a,  =  il/y/T" 


The  area  integrals  in  Equations  10  are  evaluated  similarly  using  4-point  Gaussian 
integration. 


Partial  derivatives  of  the  basis  functions  (shape  functions)  with  respect  to  the 
global  coordinates  are  evaluated  using  an  appropriate  chain  rule.  Thus, 


3<t»j  3<I>j  9xk 

9a,  3xk  9a j 
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Equations  16  can  be  solved  simultaneously  to  yield  the  partial  derivatives. 

d^j  dxm  dxn 
3<|)j  1/,2eimnerst  3ar  3as  0a, 


Expressions  for  the  partial  derivatives  with  respect  to  the  local  coordinates  appearing  in 
the  right-hand-sides  of  Equations  17  are  obtained  by  differentiating  Equations  3  and  4. 


CONSTITUTIVE  RELATIONS 

A  finite  strain  elastic-plastic  material  representation  with  von  Mises  yielding  and 
a  linear  combination  of  isotropic  and  kinematic  hardening  as  described  by  Key  and 
Krieg  is  used.2,3  The  deformation  rate  components  are 


The  normal  deformation  rate  components  are  the  material  stretchings  and  the 
mixed  components  are  the  shearings.  The  spin  components  are 


When  the  material  is  behaving  elastically,  the  stress  rate  components  are  given 
in  terms  of  the  strain  rate  and  spin  components  by 

j  a;  j  =  2/xdjj  +  X6ijdkk  +  wlkokj  +  wjkokl  (20) 

where  X  and  p  are  the  Lame  constants.  The  first  two  terms  on  the  right-hand-sides  of 
Equations  20  are  the  contributions  due  to  material  straining,  while  the  third  and 
fourth  terms  are  due  to  material  rotation. 


A  von  Mises  yield  function  is  defined  as 


0=  l$l2  -  R: 


<  0  for  elastic  behavior 
=  0  for  plastic  behavior 


(21) 


Samlia  Laboratories.  An  F.fficient  Numerical  Method  for  Time  Independent  Plasticity,  by  R.  D.  krieg. 
Albuquerque,  N.  Mcx.,  November  1977.  (SAND7 7-0943.  publication  I'NCL ASSUMED.) 
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where  R  is  the  radius  of  the  yield  surface  in  stress  space  and  l£l  is  an  effective  stress 
defined  as 


ti-y?  Mu 


(22) 


where  the  quantities  £y  are  the  deviatoric  parts  of  the  differences  between  the  stress 
components  and  the  coordinates,  <*jj,  of  the  center  of  the  yield  surface  in  stress  space; 

*-e->  £ij  =  °ij  “  ®ij  "  3^ij(°kk  "  akk^- 

The  radius  of  the  yield  surface  and  the  time  derivatives  of  the  coordinates  of 
its  center  are  given  by 


R  =  Yo+0Hep 

#  .  (^ij  —  ^ij ) 

ay  =  0  -  /3)Hep — jn - +  wikakj  +  wjkakj 


(23) 


where 


as 


Y0  =  initial  yield  stress  in  a  uniaxial  test 
E  =  elastic  modulus  from  a  uniaxial  test 
H  =  EEp/(E  -  Ep);  Ep  =  plastic  modulus  from  a  unixial  test 
)3  =  parameter  in  the  range  0  <  0  ^  1  to  control  the  hardening  behavior;  if 
0  =  0  the  behavior  is  kinematic,  if  0  =  1  the  behavior  is  isotropic. 

The  effective  plastic  strain  rate  ep  and  the  effective  plastic  strain  ep  are  defined 

•'o 


,dt 


(24) 


(The  superscript  p  denotes  the  plastic  part  of  the  deformation  rate.) 

In  the  time  integration  scheme  a  trial  stress  is  calculated  at  each  time  step 
assuming  elastic  behavior.  Thus,  using  Equations  20 


°J\  =  °!j  +  ^  Atn+ 1 12 


where  the  superscript  T  denotes  the  trial  state. 


(25) 


12 


The  trial  stress  is  substituted  into  Equation  21  to  calculate  a  trial  value  of  the 
yield  function.  If  0T  <  0  then  the  process  is  elastic  and  the  final  stress  state  is  the 
trial  stress  state.  If  0T  >  0  then  the  process  is  at  least  partly  plastic  and  further  calcu¬ 
lation  must  be  made. 


The  increment  of  effective  plastic  strain  over  the  time  step  and  the  value  of 
the  effective  plastic  strain  at  the  end  of  the  time  step  are3 


A?n+l/2 

P 


1  l£lT-Rn 
3U  1  +  H/3G 


where 


en+ 1  =  en  +  Aen+1/2 

P  p  p 

The  trial  stress  components  are  scaled  back  to  the  yield  surface  by 

0n+ 1  =aT  -  c  *T 
%  u  co«ij 

1  l£lT  -  Rn 


(26) 


(27) 


~  1  +  H/3G  |£  jT 


The  new  coordinates  of  the  center  of  the  yield  surface  are 

n+1  n  +  *n+l/2Atn+I/2 
ij  “ij  ij  “ 

where  the  time  derivatives  are  obtained  from  Equations  23  with  ?p  =  Ae”+1^/At 


(28) 

n+1/2 


CIRCUMFERENTIAL  SEGMENTATION 
REDUCTION 

Examination  of  Figure  2  reveals  that  the  method  used  to  divide  the  body  into 
elements  can  result  in  a  large  concentration  of  elements  near  the  Z-axis.  In  most 
instances,  the  region  of  primary  structural  interest,  and  the  region  where  most  elements 
are  desired,  is  the  case  wall  near  the  outer  radius.  Fine  detail  near  the  centerline  is 
superfluous  and  adds  to  the  run  time.  Moreover,  the  small  circumferential  dimensions 
of  elements  near  the  axis  can  impose  an  undesirably  small  time  step  on  the  solution. 


To  alleviate  this  problem,  a  method  for  reducing  the  number  of  elements  in 
the  circumferential  direction  by  one-half  was  devised.  This  is  illustrated  in  Figure  4. 
where  the  shaded  area  represents  the  region  of  coarser  segmentation.  At  the  transition 


FIGURE  4.  Finite  Element  Model  with  Circumferential 
Segmentation  Reduction. 
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from  finer  to  coarser  segmentation,  nodal  points  not  forming  a  corner  of  one  of  the 
coarser  elements  are  constrained  to  always  remain  midway  between  the  nodes  on  either 
side.  The  constrained  nodal  points  are  indicated  by  black  dots  in  the  figure.  The 
masses  of  the  constrained  nodal  points  are  lumped  equally  at  the  nodes  on  either  side. 
Satisfaction  of  the  virtual  displacement  equations  requires  that  the  nodal  forces  on  the 
constrained  nodal  points  also  be  divided  equally  between  the  nodes  on  either  side. 

This  technique  can  result  in  a  substantial  computational  savings  through  a 
reduction  in  the  concentration  of  elements  near  the  centerline.  It  can  also  be  used  to 
provide  a  less  detailed  representation  of  the  explosive  filler,  which  can  often  be  toler¬ 
ated  when  the  structural  integrity  of  the  case  is  of  primary  concern  in  the  analysis. 
The  suddenness  of  the  reduction  of  one-half  in  radial  segmentation  is  mitigated  by  the 
shrinking  of  the  circumferential  dimensions  of  the  elements  with  decreasing  radius  so 
that  in  practice  a  reasonably  uniform  grid  can  be  obtained.  More  than  one  level  of  seg¬ 
mentation  reduction  can  be  used  if  desired. 
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COMPUTER  IMPLEMENTATION 


THE  WHAP  CODE 

A  computer  code  called  WHAP  (Warhead  Analysis  Program)  has  been  written 
implementing  the  theory  described  in  the  previous  section.  During  problem  solution  a 
large  amount  of  nodal  and  element  data  must  be  stored  and  assessed.  It  is  not  possible 
to  fit  a  reasonable-sized  problem  into  the  standard  65K  address  space  of  the  UNIVAC 
1100/82  computer.  Two  versions  of  the  code  were  written  to  attempt  to  deal  with 
this  difficulty. 

In  the  initial  version  of  the  program  nodal  displacements,  velocities  and  forces 
were  stored  sequentially  by  plane  on  tape  (logical  unit  13).  Element  stresses,  effective 
plastic  strains,  and  yield  surface  center  coordinates  were  stored  sequentially  by  segment 
on  a  second  tape  (logical  unit  14).  Only  the  nodal  data  for  the  ith  and  (i+l)th  planes 
and  the  element  data  for  the  ith  segment  were  stored  in  core  memory  at  any  one 
time.  This  scheme  suffered  from  two  drawbacks.  First,  a  large  amount  of  input/output 
time  was  required  for  reading  and  writing  the  tapes.  Second,  it  was  not  feasible  to 
take  advantage  of  the  segmentation  reduction  scheme  described  previously. 

In  the  second  version  of  the  program,  the  “(s  FTN.O”  compiler  option  was  used 
to  permit  accessing  addresses  over  65K.  Thus,  the  entire  program  and  data  were  stored 
in  extended  memory  removing  the  need  for  storage  tapes.  Although  there  is  a  penalty 
in  terms  of  slower  execution  speed,  the  two  drawbacks  of  the  first  version  are  elim¬ 
inated.  Only  the  second  version  will  be  described  in  this  report. 


FUNCTIONING  OF  WHAP  CODE 

An  abbreviated  flow  diagram  for  the  WHAP  code  is  shown  in  Figure  5.  As  indi¬ 
cated  in  the  figure  the  functioning  of  the  program  can  be  broken  down  into  three 
phases:  setup,  solution,  and  end.  These  will  be  discussed  in  turn. 

Setup  operations  are  performed  first.  The  start  tape  for  the  problem  and  a  sol¬ 
ution  and  output  times  card  are  read.  Based  on  nodal  point  cylindrical  coordinates  in 
the  generating  plane,  the  rectangular  cartesian  coordinates  of  all  nodal  points  are 
determined.  Nodal  forces  due  to  element  surface  pressures  at  the  solution  start  arc  cal¬ 
culated.  It  is  assumed  that  the  applied  loading,  although  non-axisymmetric,  has  a  plane 
of  symmetry  such  that  only  half  the  warhead  (from  9  =  0  to  9  =  n)  need  be  anal¬ 
yzed.  The  mass  matrix  is  calculated  and  an  initial  time  step  determined.  Other  initiali¬ 
zations  as  well  as  the  printing  of  heading  information  and  solution  and  output  times 
data  are  also  carried  out. 

In  the  solution  phase,  the  approximate  time  integration  of  the  equations  of 
motion  is  performed.  At  each  time  step  the  nodal  velocities  and  positions  are  updated 
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FIGURE  S.  Flow  Diagram  of  WHAP  Code. 


using  the  central  difference  formulas,  and  new  nodal  force  contributions  due  to 
element  surface  pressures  are  calculated.  Based  on  the  new  nodal  velocities  and  posi¬ 
tions,  element  deformation  rates  and  spins  can  be  determined.  These  are  used  to 
update  the  element  stresses  through  the  constitutive  relations.  Nodal  force  contribu¬ 
tions  due  to  element  stresses  are  calculated  and  added  to  the  forces  due  to  surface 
pressures.  With  the  integration  step  essentially  complete,  an  output  data  print  and/or 
plot  tape  writes  can  be  made  if  called  for  by  the  output  times  card.  Finally,  checks 
are  made  to  see  if  the  end  of  the  solution  time  period  has  been  reached  or  if  an 
excessive  time  step  drop  has  occurred  indicating  that  the  solution  should  be  termi¬ 
nated.  If  not,  the  integration  procedure  is  repeated. 

At  the  end  of  the  run  a  restart  tape  having  the  same  format  as  the  original 
start  tape  is  written  (except  that  additional  nodal  and  element  data  are  included). 
Problem  solution  can  be  continued  using  this  tape  and  a  new  solution  output  times 
card.  This  feature  permits  checking  to  see  how  a  problem  is  progressing  before  too 
much  computer  time  has  been  wasted  on  an  unsatisfactory  solution. 


OVERALL  PROBLEM  SOLUTION 

The  overall  problem  solution  procedure  including  the  pre-  and  post-processing 
required  for  use  of  the  WHAP  code  is  indicated  in  the  flow  diagram  in  Figure  6.  A 
two-dimensional  quadrilateral  mesh  for  the  generating  plane  is  created  interactively  on 
the  HP9845A  Desktop  Computer  using  program  MESH.  Interactive  input  of  other 
required  data  is  accomplished  using  program  INWHAP.  These  mesh  and  input  data  are 
transmitted  in  the  form  of  card  images  to  the  UNIVAC  1100/82  computer  and  stored 
in  a  File  using  the  HP9845A  as  a  remote  terminal.  Program  WHAPST  is  then  executed 
on  the  UNIVAC  1100/82  to  transform  the  file  data  into  a  start  tape  for  the  WHAP 
program.  A  small  amount  of  additional  data  (regarding  solution  and  output  times)  is 
input  during  running  of  the  WHAP  code  itself,  which  performs  the  approximate  time 
integration.  Printed  output  and  plot  tapes  are  generated  during  execution.  Post-process¬ 
ing  programs  DFSTRW  and  CNTRW  are  used  to  obtain  deformed  structure  and  stress 
contour  plots  from  the  plot  tapes.  Additional  plot  specification  cards  are  required  for 
the  running  of  these  programs. 


USER  INFORMATION 

For  reference  purposes,  a  list  of  most  of  the  variables  used  in  the  WHAP  code 
is  given  in  Appendix  A.  Brief  descriptions  of  the  subroutines  making  up  the  code  are 
given  in  Appendix  B.  User  instructions  for  the  WHAP  code  and  associated  pre-  and 
post-processors  are  given  in  Appendix  C. 


Times 
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FIGURE  6.  Flow  Diagram  of  Overall  Problem  Solution. 
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CODE  CHECKOUT 


ways: 


Checkout  of  a  computer  code  such  as  this  can  be  accomplished  in  two  main 


1.  By  using  the  code  to  solve  a  problem  for  which  an  analytic  solution  exists 

2.  By  comparing  the  code  results  with  results  of  another,  already-verified  code 
for  the  same  problem. 


The  first  of  these  methods  is  doubtless  the  more  elegant.  However,  in  the  case 
of  bodies  subjected  to  non-axisymmetric  impulsive  loadings,  suitable  analytic  solutions 
are  not  available.  Hence,  the  second  method  was  used  (although  an  attempt  was  made 
to  indicate  that  the  results  are  reasonable  in  the  light  of  elementary  theory).  The 
HONDO  11  code  was  used  for  comparison  purposes.  Two  checkout  problems  were  run. 


Checkout  Problem  No.  1 

The  first  checkout  problem  involved  an  elastic-plastic  rod  subjected  to  a  step 
pressure  loading  at  one  end.  The  rather  coarse  finite  element  grid  used  for  both  the 
WHAP  and  HONDO  II  runs  is  shown  in  Figure  7.  The  geometric  and  material  param¬ 
eters  defining  the  problem  are  also  indicated.  Eight  circumferential  segments  were  used 
in  the  WHAP  code  model. 


mum 
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L/D  =  5 
E/Ep  -  4 
p  s  .3 
p/Y  s  2 


NWC  TP  6369 


Axial  stress  versus  distance  from  the  impacted  end  of  the  rod  is  plotted  in 
Figure  8.  The  HONDO  II  and  WHAP  results  are  quite  close  to  one  another.  The  use 
of  more  segments  would  probably  have  produced  even  better  agreement.  The  dashed 
line  in  the  figure  is  the  result  obtained  from  one-dimensional  rod  theory.  Differences 
between  the  finite  element  and  theoretical  results  are  due  primarily  to  transverse 
inertia  effects.  A  rod  with  a  higher  length-to-diameter  ratio  might  have  yielded  results 
more  nearly  in  accord  with  one-dimensional  theory.  However,  errors  in  the  code  associ¬ 
ated  with  transverse  effects  might  not  have  been  apparent  in  such  a  model. 


DISTANCE  FROM  FRONT  END,  z/L 

FIGURE  8.  Axial  Stress  Versus  Distance  From  Front 
End  for  Elastic-Plastic  Rod  Subjected  to  Step 
Pressure. 


Checkout  Problem  No.  2 

The  loading  in  the  first  checkout  problem  was  axisymmetric.  Hence,  this  first 
problem  did  not  test  the  non-axisymmetric  part  of  the  coding.  The  second  checkout 
problem  involved  an  elastic  ring  subjected  to  a  ramp  pressure  on  its  outer  surface  that 
varied  sinusoidally  around  the  circumference  of  the  ring.  Plane  strain  was  assumed. 
The  HONDO  II  code  has  the  capability  of  solving  such  a  plane  strain  problem.  (A 
large  constant  value  is  added  to  all  the  nodal  radii  in  the  model.)  A  WHAP  code  solu¬ 
tion  was  obtained  treating  the  ring  as  a  solid  of  revolution  subjected  to  a  non-axisym- 
metric  loading.  The  finite  element  models  used  and  the  parameters  defining  the 
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problem  are  given  in  Figure  9.  One  hundred  circumferential  segments  were  used  in  the 
WHAP  model,  which  is  the  same  as  the  number  of  quadrilaterals  around  the  circum¬ 
ference  of  the  HONDO  II  model. 

Stress  parallel  to  the  surface  in  the  first  layer  of  elements  beneath  the  load  is 
plotted  versus  circumferential  angle  in  Figure  10.  The  slight  difference  between  the 
WHAP  and  HONDO  II  results  may  be  due  to  a  difference  in  the  way  the  artificial 
viscosity  is  calculated  in  the  two  codes.  In  an  axisymmetric  elastic  body  subjected  to 
sinusoidal  loading  one  would  expect  sinusoidal  variation  in  the  stress  as  well.  Except 
for  the  slight  bump  at  low  angles,  the  curves  are  nearly  sinusoidal.  (A  sine  curve  with 
the  same  maximum  amplitude  could  have  been  superposed  on  the  finite  element 
results  but  was  omitted  for  clarity.) 


HONDO  II  grid  WHAP  grid 


FIGURE  9.  Elastic  Ring  Subjected  to  Ramp  Pressure 
Sinusoidally  Varying  Around  Circumference. 
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CIRCUMFEPENTIAL  ANGLE,  0.  deg 

FIGURE  10.  Circumferential  Stress  vs.  Angle  for  Elastic  Ring 
Subjected  to  Ramp  Pressure. 


ANALYSIS  OF  PROJECTILE  IMPACTIi  G  SIMULATED  CONCRETE 


SMALL  PROJECTILE  FIRINGS 

Experimental  firings  of  small  steel  projectiles  at  normal  incidence  against  sim¬ 
ulated  concrete  half-space  and  thin  steel  plate  targets  are  described  by  Stronge  and 
Schulz.4  The  uninstrumented  test  projectiles  were  fiat-fronted  cylinders  with  an  inter¬ 
nal  cavity  and  were  intended  to  roughly  represent  warheads.  During  penetration  the 
projectiles  developed  a  bulge  near  the  front  of  the  internal  cavity  and.  at  high  enough 
impact  velocities,  fractured  in  this  region.  Projectile  deformations  measured  after  testing 
were  compared  with  finite  element  analysis  results  obtained  using  the  HONDO  II  code. 

Similar  firings  of  projectiles  against  targets  at  obliquity  are  described  by  Schulz 
and  Heimdahl.5  Some  of  the  projectiles  tested  at  obliquity  were  filled  with  an  explo¬ 
sive  simulant  (plaster  of  Paris)  while  others  were  left  unfilled.  A  cross-sectional  view  of 


4  W.  J.  Stronge  and  J.  C.  Schulz.  “Projectile  Impact  Damage  Analysis,"  in  Papers  Presented  at  the  Syntpo- 
sium  on  Computational  Methods  in  Son-linear  Structural  and  Solid  Mechanics ,  edited  by  Almicd  K.  Noor  and 
Harvey  G.  McComb,  Jr.  Pcrgamon  Press,  New  York,  1981.  (Also  published  as  special  issue  of  Journal  of  Computers 
and  Structures,  Vol.  13,  No.  1-2  (2981).  op.  287-294.) 

®  Naval  Weapons  Center.  Oblique  Impact  Projectile  Damage,  by  J.  C.  Schulz  and  O.  I  .  R,  Heimdahl. 
China  l.ake,  Calif.,  NWC,  January  1982.  (NWC  TM  4695,  publication  UNCLASSIFIED.) 
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the  projectile  used  in  oblique  impact  firings  is  shown  in  Figure  11.  In  addition  to 
bulging,  projectiles  fired  at  obliquity  also  bent  such  that  the  centerline  of  the  nose  and 
case  portions  were  not  colinear. 

Measurements  of  the  final  deformed  shapes  of  projectiles  from  these  oblique 
impact  firings  were  compared  with  WHAP  code  analytical  results.  Such  a  comparison 
does  not  provide  an  especially  critical  test  of  code  performance.  Even  when  reasonable 
agreement  is  obtained,  stress  values  may  differ  considerably  between  analysis  and 
experiment.  However,  since  projectile  stresses  were  not  measured  (due  to  the  difficul¬ 
ties  in  providing  suitable  instrumentation),  no  other  comparison  appeared  feasible. 


PROJECTILE  ANALYSIS  TECHNIQUE 

A  structural  analysis  technique  for  projectiles  impacting  hard  targets  at  normal 
incidence  has  been  developed  at  the  Naval  Weapons  Center.1,4  This  technique  involves 
finding  an  approximate  loading  history  to  be  applied  to  the  front  of  the  projectile  to 
represent  forces  at  the  target-projectile  interface  and  then  using  this  approximate 
loading  in  a  nonlinear  finite  element  analysis  to  calculate  stresses  and  deformations  in 
the  projectile.  This  same  technique  can  be  used  for  projectiles  impacting  at  obliquity  if 
the  loading  history  is  modified  to  account  for  non-axisymmetric  effects. 
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PROJECTILE  SELECTED  FOR  ANALYSIS 

To  keep  run  costs  down,  an  unfilled  projectile  (without  the  plaster  of  Paris) 
was  considered  for  analysis.  The  impact  velocity  of  the  projectile  was  taken  as  2,155 
ft/s  (corresponding  to  one  of  the  unfilled  test  projectiles  fired  at  normal  incidence). 
WHAP  code  runs  were  made  using  loading  histories  corresponding  to  penetration  into  a 
concrete  half-space  at  normal  incidence  and  at  45  degrees  obliquity.  For  comparison 
purposes,  a  HONDO  II  run  at  normal  incidence  was  also  made.  The  solution  was 
stopped  at  a  maximum  of  75  microseconds  by  which  time  plastic  deformation  of  the 
projectile  appeared  to  have  ceased. 


FINITE  ELEMENT  MODEL 

A  model  of  the  projectile  containing  264  nodal  points  and  198  elements  in  the 
generating  plane  was  constructed.  Ten  radial  planes  were  used  to  divide  the  projectile 
into  nine  segments  circumferentially.  The  4340  steel  was  modeled  as  an  elastic, 
perfectly  plastic  material  with  the  properties  given  in  Table  1. 


TABLE  1.  Material  Properties  for  4340  Steel. 


Young's  modulus,  psi .  30,000.000 

Poisson's  ratio  .  0.3 

Dynamic  yield  stress,  psi .  177,000 

Density,  lb  sec2/in4  .  0.000733 


LOADING  HISTORIES 

For  impact  at  normal  incidence  a  loading  history  was  obtained  using  the 
method  of  Stronge  and  Schulz.4  This  loading  history  consisted  of  a  quasi-steady-state 
penetration  theory  drag  force6  to  which  was  added  an  initial  transient  to  account  for 
start-up  conditions.  The  force  was  applied  as  a  pressure  uniformly  distributed  over  the 
front  end  of  the  projectile.  For  impact  at  obliquity  two  different  loading  histories 
were  used. 

In  the  first  loading  (Load  No.  I),  the  same  pressure  as  a  function  of  time  was 
used  as  for  normal  impact,  except  that  the  application  time  was  shifted.  Initially, 
pressure  was  applied  only  at  the  leading  edge  of  the  projectile  (the  edge  that  impacted 
the  target  first).  As  the  projectile  moved  forward  and  more  of  its  front  end  crossed 
the  plane  of  the  target  the  pressure  was  applied  over  this  portion  of  the  front  end  as 


**  Defense  Nuclear  Agency.  Projectile  Penetration  in  Harth  Materials  Theory  and  Computer  Analysis,  be  R. 
S.  Bernard  and  D.  (\  Creighton.  Washington,  D.C.,  November  1977.  (Contract  Report  S-76-13.  puhli  ti on 
UNCLASSIFIED.) 
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well.  Only  when  the  projectile  was  completely  imbedded  was  pressure  applied  over  the 
entire  front  end.  The  effect  was  a  pressure  that  swept  across  the  front  surface  of  the 
projectile  as  this  surface  reached  the  target. 

This  loading  did  not  produce  results  in  agreement  with  experiment.  Consequent¬ 
ly,  a  second  loading  (Load  No.  2)  was  used.  In  this  second  loading,  the  magnitude  of 
the  pressure  at  the  leading  edge  was  increased  by  5rr  and  the  magnitude  of  the 
pressure  at  the  trailing  edge  was  decreased  by  5%.  The  pressure  was  allowed  to  vary 
linearly  across  the  front  surface  between  these  extremes.  Some  rationale  for  this 
second  loading  will  be  given  later. 


FINITE  ELEMENT  RESULTS 

The  WHAP  code  results  for  the  final  deformed  shapes  of  projectiles  impacting 
simulated  concrete  at  normal  incidence  and  at  45  degrees  obliquity  (Load  No.  2)  are 
shown  in  Figure  12.  Qualitatively,  the  analytically  determined  shapes  agree  well  with 


FIGURE  12.  Final  Deformed  Shapes  of  WHAP  Models  of  Projectiles  Impacting  Simulated 
Concrete  at  2,155  ft/s  at  Normal  Incidence  and  45  Degrees  Obliquity. 
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the  shapes  of  test  projectiles  illustrated  photographically  by  Schulz  and  Heimduhl.5  To 
obtain  a  quantitative  comparison,  the  increase  in  diameter  at  the  bulge,  the  angle 
between  nose  and  case  centerlines  and  the  distance  from  the  bulge  to  the  rear  of  the 
projectile  (D,  d,  and  L  in  the  figure)  were  measured  for  the  test  projectiles.  A  tabula¬ 
tion  of  measured  and  analytically  determined  values  is  given  in  Table  2.  (No  firing  was 
made  at  2,155  ft/s  and  45  degrees  obliquity  so  data  at  2,200  and  2.240  ft/s  were 
3  extrapolated  for  the  table.) 


TABLE  2.  Comparison  of  Analytical  and  Test  Results  for  Small  Projectile  Impacting 
Simulated  Concrete  at  2,155  ft/s. 


1 

Normal  incidence 

45  degrees  obliquity 

HONDO  II 

WHAP 

Test 

projectile 

WHAP 

Load  no.  1 

WHAP 

Load  no.  2 

Test 

projectile 

Increase  in  diameter  at 

bulge,  AO,  in. 

0.088 

0.046 

0.047 

0  052 

Angle  between  nose  and 

case  centerlines.  0.  deg. 

0.9 

3.7 

4.4 

Distance  from  bulge  to 

rear  of  projectile,  L,  in. 

1.60 

1.59 

1.59 

For  impact  at  normal  incidence  there  . ;  good  agreement  between  the  predic¬ 
tions  of  the  HONDO  II  and  WHAP  codes  and  t.  e  measured  values.  For  impact  at  45 
degrees  obliquity,  the  diameter  increases  and  the  distances  from  the  bulge  to  the  rear 
of  the  projectiles  predicted  by  the  WHAP  code  with  the  two  different  loadings  are  in 
good  agreement  with  each  other  and  also  with  the  experimental  results.  However,  the 
I  bend  angle  predicted  with  Load  No.  1  is  considerably  less  than  that  predicted  with 

Load  No.  2,  which  agrees  fairly  well  with  experiment.  An  explanation  for  this  behavior 
is  as  follows:  Compared  to  impact  at  normal  incidence,  there  is  extra  target  material 
on  the  leading  edge  side  of  the  projectile.  This  extra  material  can  be  expected  to  pro¬ 
duce  some  additional  confinement  and  an  increase  in  the  pressure  loading  on  this  side. 
I  Conversely,  there  is  a  lack  of  material  on  the  trailing  edge  side  of  the  projectile,  which 

can  be  expected  to  reduce  the  pressure  loading.  The  ±5%  adjustment  that  was  made  to 
the  loading  magnitude  has  no  especial  justification  except  that  it  produced  better 
results. 


RUN  COSTS 

Oblique  impact  runs  for  this  problem  using  the  WHAP  code  proved  to  be  quite 
expensive.  Run  costs  on  the  Naval  Weapons  Center  UNIVAC  1  100/82  were  approx¬ 
imately  $2,700  versus  $140  for  comparable  WHAP  and  HONDO  II  runs.  CPU  times 
per  time  step  were  6.67  versus  0.34  with  the  two  programs. 
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CONCLUSION 

A  program  for  the  structural  analysis  of  axisymmetric  warheads  impacting 
targets  at  obliquity  has  been  developed.  The  WHAP  code  itself,  as  well  as  pre-  and 
post-processors  required  for  complete  problem  solution,  have  been  written  and  checked 
out,  and  the  entire  package  is  now  operational.  Results  obtained  for  the  analysis  of  a 
small  steel  projectile  penetrating  into  a  simulated  concrete  half-space  at  both  normal 
incidence  and  45  degrees  obliquity  appear  reasonable. 

One  objective  of  this  work  was  to  produce  a  code  that  is  inexpensive  to  use. 
This  goal  was  not  attained.  Run  costs  using  the  WHAP  code  are  an  order  of  magnitude 
higher  than  for  the  two-dimensional  code,  HONDO  II.  As  a  result,  analyses  of  war¬ 
heads  impacting  at  obliquity  using  the  new  code  will  probably  not  be  performed  as 
routinely  as  analyses  of  warheads  impacting  at  normal  incidence.  However,  the  code  is 
available  for  situations  where  the  effects  of  obliquity  are  of  especial  interest  -and 
where  money  is  available  for  making  runs. 

The  small  projectile  problem  that  was  run  points  out  the  need  for  better  char¬ 
acterization  of  the  loading  histories  of  projectiles  impacting  targets  at  obliquity. 
Reasonable  agreement  between  code  and  test  results  was  obtained  only  by  artificially 
adjusting  the  loading  used.  Almost  no  data  regarding  the  effects  of  obliquity  on  pene¬ 
tration  forces  are  available,  especially  for  the  initial  transient  stage  of  loading  behavior. 
Work  needs  to  be  done  in  this  area  to  provide  input  for  running  of  the  WHAP  code. 
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Appendix  A 

DEFINITIONS  OF  VARIABLES 

Definitions  of  most  of  the  variables  used  in  the  WHAP  program  are  listed  for 
reference  in  this  appendix. 

COMMON/ZVAR/variables 

Z(NDIMZ),  main  storage  array  for  program 
IZ(NDIMZ),  integer  equivalent  of  Z  array 
NDIMZ,  maximum  dimension  of  Z  array 

N2-N22,  integers  giving  starting  locations  for  program  variables  stored  in  Z  array 
COMMON/BK1 /variables 

NNP,  number  of  nodal  points  in  generating  plane 
NEL,  number  of  quadrilateral  element  traces  in  generating  plane 
MPLN,  maximum  number  of  radial  planes  in  circumferential  direction  (in  regions 
of  finest  circumferential  segmentation) 

MSEG,  maximum  number  of  elements  in  circumferential  direction  =  MPLN  1 
MNP,  total  number  of  nodal  points 
MEL,  total  number  of  elements 

NCNP,  number  of  circumferential  arcs  with  constrained  nodal  points  (for  transi¬ 
tion  to  coarser  circumferential  segmentation) 

NMAT,  number  of  materials 

NESP,  number  of  element  surface  pressures 

NLTC,  number  of  load-time  curves 

MLP,  maximum  number  of  load  points  on  any  load-time  curve 

HED(20),  problem  title 

IDT,  date  initial  start  tape  created 

ITM,  time  of  day  initial  start  tape  created 

NSTRT,  start  number 

VZO,  initial  axial  velocity 

COMMON/BK2/variablcs 

A(8),  values  of  local  a-coordinate  at  Gaussian  quadrature  points 

B(8),  values  of  local  b-eoordinate  at  Gaussian  quadrature  points 

C(8),  values  of  local  c-coordinate  at  Gaussian  quadrature  points 

P(4,4),  shape  function  values  at  Gaussian  quadrature  points  in  element  surface 
used  in  calculation  of  loads  due  to  surface  pressures 
PP(8,8),  shape  function  values  at  Gaussian  quadrature  points  inside  element 
DPDA(8,8),  derivatives  of  shape  functions  with  respect  to  a-coordinate  at  Gaussian 
quadrature  points 

DPDB(8,8).  derivatives  of  shape  functions  with  respect  to  b-coordinate  at  Gaussian 
quadrature  points 

DPDC(8,8),  derivatives  of  shape  functions  with  respect  to  c-coordinate  at  Gaussian 
quadrature  points 

THSEG,  tr/MSEG,  angle  between  planes  corresponding  to  finest  segmentation 
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COM  M  ON/BK3/variables 

TT,  current  time  in  solution  integration 
DTI,  old  time  step 
DT2,  new  time  step 

DTO,  initial  time  step  x  10-5,  used  to  check  for  time  step  drop 

DTN,  time  step  in  a  given  element,  used  in  determining  new  time  step 

NELCRT,  element  with  critical  time  step 

NSGCRT,  segment  containing  element  with  critical  time  step 

NSTEP,  number  of  current  integration  step 

Bl,  coefficient  of  quadratic  term  in  artificial  viscosity 

B2,  coefficient  of  linear  term  in  artificial  viscosity 

B3,  coefficient  in  keystone  viscosity 

COMMON/BK4/variables 

TPRT1,  initial  time  for  printed  output 

TPRT2,  final  time  for  printed  output,  solution  end 

DTPRT,  time  interval  for  printed  output 

TPLT1,  initial  time  for  plotted  output 

TPLT2,  final  time  for  plotted  output 

DTPLT,  time  interval  for  plotted  output 

NDOUT(50),  nodes  at  which  displacement  quantities  are  to  be  output 
NLOUT(50),  elements  in  which  stresses  are  to  be  output 
NPLOUT(20),  planes  for  which  displacement  quantities  are  to  be  output 
NSGOUT(20),  segments  for  which  stresses  ax  to  be  output 

COMMON/BK5/variables 

DF,  determinant  of  deformation  gradient 
D(6),  deformation  (strain)  rate 
W(3)  spin  (non-zero  components) 

Quantities  stored  unchanged  in  Z  array  throughout  solution 

Z(1)~R(NNP),  radial  coordinates  of  nodes  in  undeformed  body 
Z(N2)~Z(NNP),  axial  coordinates  of  nodes  in  undeformed  body 
IZ(N3)~BC(NNP),  nodal  boundary  condition  codes 

IZ(N4)~NPLN(NNP),  number  of  radial  planes  associated  with  nodal  points  in 
generating  plane 

IZ(N5P"ND(NNP),  global  nodal  numbers  cooresponding  to  generating  plane  nodal 
numbers 

IZ(N6)~NDA(NNP),  constrained  node  indicator  for  use  in  displacement  update 
lZ(N7pNDC(NCNP),  nodal  numbers  associated  with  circumferential  ares  contain¬ 
ing  constrained  nodes 

1Z(N8)~NSEG(NEL),  number  of  segments  associated  with  element  traces  in 
generating  plane 

1Z(N9)~NL(NEL),  global  element  numbers  corresponding  to  generating  plane 
element  numbers 

IZ(N10)^IX(5,NEL),  element  nodal  point  and  material  specifications 
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Z(N  1  I  )~PR0P(15,NMAT),  material  properties 

IZ(N12)~ISP(5,NESP),  nodal  points  and  load-time  curves  for  pressurized  surfaces 
Z(N13)~PM(4,NESP),  pressure  multipliers  at  nodes  on  pressurized  surfaces 
Z(N  14)~AT(NESP),  arrival  times  for  surface  pressures 
Z(N15)~LT(2,NLTC,MLP).  load-time  points 
Z(N16)~MS(NNP),  nodal  masses 

Z(NI7)~DJ0(NEL),  initial  Jacobian  determinants  at  element  centers 

Nodal  quantities  updated  in  Z  array  throughout  solution 

Z(N18)~XX(3,MNP),  coordinates  of  nodal  points  in  deformed  body 
Z(NI9)~VV(3,MNP),  velocity  components  of  nodal  points  in  deformed  body 
Z(N20)~FF(3,MNP),  force  components  at  nodal  points  in  deformed  body 

Stress  quantities  updated  in  Z  array  throughout  solution 
Z(N21)~S(6,MEL),  element  stress  components 

Z(N22)~STRG(6,MEL),  storage  for  element  yield  surface  centers,  effective  plastic 
strains  and  state  variables 

Material  properties  (1—6  are  input,  7—15  are  calculated) 

PROP(l,NMAT)~RO,  density 
PROP(2,NMAT)~E,  Young’s  modulus 
PROP(3,NMAT)''NU,  Poisson’s  ratio 
PROP(4,NMAT)~SY,  uniaxial  yield  stress 
PROP(5,NMATrET,  hardening  modulus 
PROP(6,NMAT)~"BETA,  hardening  parameter 
PROP(7,NMAT)~LMBDA=TWOG*NU/(l.-2.*NU) 

PROP(8,NMAT)~TWOG=E/(l  ,+NU) 

PROP(9,NMAT)~H=E*ET/(E-ET) 

PROP(  1 0,NM  AT)~YR=. 8 1 649658 1  *SY 
PROP(  1 1  ,NMAT)~CYR=.8 1 649658 1  *BETA*H 
PROP(  1 2,NMAT)~CEP=.8 1 6496581  *TWOG 
PROP(  1 3,NM  AT)~CCIG=  1 ./( 1  +H/ 1 ,5/TWOG) 
PROP(14,NMAT)~CCA=H/1.5/TWOG*(l.-BETA) 
PROP(15,NMAT)~CX=SQRT((LMBDA+TWOG)/RO) 


31 


NWC  TP  6369 


Appendix  B 

SUBROUTINE  DESCRIPTIONS 

Brief  descriptions  of  the  subroutines  making  up  the  WHAP  code  are  given  in 
this  appendix. 

MAIN  PROGRAM  WHAP 

This  short  routine  acts  as  a  driver  for  the  remainder  of  the  program.  The  size 
of  the  common  variable  Z,  where  most  data  is  stored,  can  be  changed  by  changing  the 
parameter  NDIMZ  in  this  routine.  WHAP  calls  subroutines  SETUP  and  SOLVE. 

SUBROUTINE  SETUP  (N2,N3,N4,N5,N6,N7,N8,N9,N10,N1 1  ,N1  2,Nl 3,N14, 

N 1 5  ,N  1 6,N  1 7,N  1 8,N  1 9,N20,N2 1  ,N22,  NDIMZ) 

Reading  of  the  start  tape,  and  storage  and  initialization  of  problem  variables  are 
accomplished  in  this  subroutine.  SETUP  calls  subroutines  LOAD,  MASS,  DETJ  and 
STEP. 

SUBROUTINE  LOAD  (ISP,PM,AT,LT,XX,FF,TT) 

Nodal  forces  due  to  element  surface  pressure  loads  are  calculated  in  this  sub¬ 
routine.  Gaussian  quadrature  is  used  to  evaluate  the  area  integrals  involved. 

SUBROUTINE  MASS  (XX,ND,NDA,NSEG,IX,RO,MS) 

Nodal  masses  are  calculated  in  this  subroutine.  Gaussian  quadrature  is  used  to 
evaluate  the  volume  integrals  involved.  An  adjustment  in  the  case  of  nodes  at  the  tran¬ 
sition  from  finer  to  coarser  segmentation  is  performed. 

SUBROUTINE  DETJ  (XX.IX.DJ) 

The  initial  Jacobian  determinant  of  the  transformation  of  Equations  4  is  calcu¬ 
lated  at  the  center  of  each  element. 

SUBROUTINE  STEP  (XX,1X,PR0P) 

An  initial  time  step  that  ensures  solution  stability  is  calculated.  At  subsequent 
times  this  calculation  is  performed  in  subroutine  FORCE. 

SUBROUTINE  SOLVE  fR,Z,BC'.NPLN,ND,NDA,NDC,NSEG,NL,lX,PROP, 

ISP.PM.AT.LT,  MS,DJO,XX.VV,FF,S.STRG) 

The  approximate  time  integration  solution  is  carried  out  in  this  subroutine. 
Nodal  velocities  and  positions  are  updated  and  new  nodal  forces  due  to  surface 
pressure  loads  are  calculated.  In  a  loop  over  the  number  of  elements,  element  deforma¬ 
tion  rates  and  spins  are  calculated  and  used  to  update  the  element  stresses.  Contribu¬ 
tions  to  the  nodal  forces  due  to  element  stresses  are  calculated  and  a  new  solution 
time  step  is  determined.  An  output  data  print  and  plot  tape  writes  are  accomplished 
when  requested.  After  the  solution  is  completed  a  restart  tape  is  written.  SOLVE  calls 
subroutines  LOAD,  STRAIN,  STRESS  and  FORCE. 
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SUBROUTINE  STRAIN  (ND1  ,ND2,ND3,ND4,ND5,ND6.ND7.ND8,XX,VV.DJ0) 

Calculation  of  the  deformation  rate  and  spin  at  the  element  center  is  accom¬ 
plished  in  this  subroutine.  The  deformation  gradient  at  the  element  center  is  also 
determined. 

SUBROUTINE  STRESS  (S.LMBDA,TWOG.YR,CYR,CHP.CCSIG.CCA,A,EP,STATE) 

The  constant  element  stress  based  on  the  deformation  rate  at  the  element 
center  is  calculated.  An  elastic  trial  stress  is  calculated  and  then  scaled  back  to  the 
yield  radius  if  plastic  flow  has  occurred. 

SUBROUTINE  FORCE  (ND1,ND2(ND3>ND4)ND5,ND6,ND7,ND8> 

RO,CX,XX,VV.FF,S,ISEG,NSEG) 

Nodal  forces  due  to  the  element  stress  are  determined.  Artificial  and  keystone 
viscosity  contributions  are  included.  An  element  time  step  is  calculated. 
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Appendix  C 

USER  INSTRUCTIONS 

User  instructions  for  the  input  of  data  to  the  pre-processor  WHAPST,  to  the 
WHAP  program  itself,  and  to  the  post-processors  DFSTRW  and  CNTRW  are  provided 
in  this  appendix. 

A.  Pre-processor  WHAPST 

This  program  generates  a  start  tape  (logical  unit  8)  for  the  WHAP  code  and 
also  creates  a  plot  of  the  undeformed  generating  plane.  Input  data  for  this  program  is 
usually  created  interactively  on  the  HP9845A  Desktop  Computer  and  transmitted  in 
the  form  of  card  images  to  the  UNIVAC  1 100/82.  However,  these  card  images  could 
also  be  obtained  by  other  means.  The  format  for  these  card  images  is  given  below. 

1.  Heading  card  (20A4) 

Columns  1-80  HED(20),  information  to  identify  run 

2.  Mesh  file  card  (A6) 

Columns  1—6  MF,  name  of  mesh  file  on  HP9845A  (can  be  left  blank 
if  data  not  generated  on  HP9845A) 

3.  Storage  file  card  (A6) 

Columns  1  -6  SF,  name  of  storage  file  on  HP9845A  (can  be  left 
blank) 


Control  card 
Columns  1-: 
Columns  6- 
Columns  1 1 
Columns  16- 
Columns  21- 
Columns  26- 
Columns  31- 
Columns  36- 

Columns  41  - 
Columns  46- 
Columns  5 1  - 


(1015, El  0.0) 


NNP,  number  of  nodal  points  in  generating  plane 

NEL,  number  of  elements  in  generating  plane 

MPLN,  maximum  number  of  radial  planes 

MSEG  =  MPLN— 1 

NMAT,  number  of  materials 

NESP,  number  of  element  surface  pressures 

NLTC,  number  of  load-time  curves 

MLP,  maximum  number  of  points  on  any  load-time 

curve 

NNC,  number  of  nodal  point  specification  cards 
NEC,  number  of  clement  block  specification  cards 
VZ0,  initial  axial  velocity 


5.  Nodal  output  cards  (1015,  5  cards) 

List  of  nodal  points  for  which  printed  output  is  desired  (50  points 
maximum,  first  zero  or  blank  terminates  list) 

6.  Element  output  cards  (1015,  5  cards) 

List  of  elements  for  which  printed  output  is  desired  (50  elements 
maximum) 
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7.  Radial  plane  output  cards  (1015,  2  cards) 

List  of  radial  planes  for  which  printed  output  is  desired  (20  planes 
maximum) 

8.  Segment  output  cards  (1015,  2  cards) 

List  of  segments  for  which  printed  output  is  desired  (20  segments 
maximum) 


9. 


Nodal  point 
Columns  1- 
Columns  6- 
Columns  16 
Columns  26 
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Columns  35 
Columns  41 
Columns  51 


specification  cards  (15,2E10.0,15,4X,13.3X.2L10.0,  NNC  cards) 

5  N,  nodal  point  number 

15  R(N),  radial  coordinate 

Z(N),  axial  coordinate 

INCR,  nodal  point  generation  increment,  taken  from  last 
card  read,  BC  of  intermediate  points  is  from  previous 
card 

BC(N),  displacement  boundary  condition  code 
SR,  spacing  ratio  for  incremented  points 
RAD,  radius  of  arc  connecting  end  points  along  which 
points  are  generated 
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10.  Element  block  specification  cards  (915,  NEC  cards) 


Columns  1-5 

Columns  6-10 

Columns  11-15 
Columns  16-20 
Columns  21-25 
Columns  26-30 

Columns  31—35 
Columns  36-40 
Columns  41-45 


MCI,  nodal  point  number  of  one  comer  node  of 
four-sided  element  block 

MC2,  nodal  pom:  number  of  second  corner  node  of 
four-sided  element  block  (progressing  counter-clockwise) 
MC3,  nodal  point  number  of  third  corner  node 
MC4,  nodal  point  number  of  fourth  comer  node 
MATL,  material  identification  number 
NLGRP,  element  group  number  for  determining  number 
of  segments,  if  NLGRP  =  0  then  NLGRP  =  1 
NELRW,  number  of  elements  in  row 
NROWS,  number  of  rows 

NGEN,  generation  parameter,  if  NGEN  =  0  nodal  points 
in  interior  of  block  are  found  as  intersections  of  straight 
lines  connecting  nodal  points  on  corresponding  sides,  if 
NGEN  >  0,  interior  nodal  points  are  found  using 
Laplace  generation  with  NGEN  iterations 


11.  Material  specification  cards  (NMATL  sets) 

a.  Title  card  (20A4) 

Columns  I  80  MT1TLE,  title  to  identify  material 

b.  Property  card  (6E10.0) 

Columns  1-10  PROP(l),  material  density,  RO 
Columns  11-20  PROP(2),  Young’s  modulus,  E 
Columns  21-30  PROP(3),  Poisson’s  ratio,  NU 
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Columns  3 1  -40  PROP(4),  yield  stress.  SY 
Columns  41-50  PROP(5),  hardening  modulus,  ET 
Columns  51-60  PROP(6),  hardening  parameter,  BETA 

12.  Element  surface  pressure  specification  cards  (415,  SE10.0,  NESP  cards) 
Columns  1-5  ISP(l),  first  nodal  point 

Columns  6—10  ISP(2),  second  nodal  point,  two  points  must  form  one 

edge  of  element  in  generating  plane 

Columns  11  —  15  ISP(3),  number  of  segment  in  circumferential  direction 

Columns  16-20  ISP(4),  number  of  load-time  curve 

Columns  21-30  PM(1),  pressure  multiplier  at  first  corner  of  surface 

element 

Columns  31-40  PM(2),  pressure  multiplier  at  second  comer 

Columns  41-50  PM(3),  pressure  multiplier  at  third  corner 

Columns  51-60  PM(4),  pressure  multiplier  at  fourth  comer 

Columns  61-70  AT,  arrival  time  of  pressure  at  element  surface 

13.  Load-time  curve  specification  cards  (NLTC  sets) 
a.  Load-time  cards  (2E10.0,  MLP  cards) 

Columns  1-10  LT(1),  time 

Columns  11-20  LT(2),  load 

B.  WHAP  program 

Most  input  data  for  this  program  is  obtained  from  the  start  tape  written  by 

WHAPST  (or  from  a  restart  tape  written  by  WHAP  itself).  This  start  tape  is  read  by 

WHAP  on  logical  unit  8.  Three  blank  tapes  must  be  assigned  for  the  writing  of  a 

restart  tape,  a  displacement  data  plot  tape  and  a  stress  data  plot  tape  on  logical  units 

9,  20,  and  25,  respectively.  A  small  amount  of  additional  data  regarding  solution  and 
output  times  must  be  input  in  the  format  given  below. 

1.  Solution  and  output  times  card  (6F10.0) 

Columns  1  — 10  TPRT1,  first  time  for  printed  output 
Columns  11-20  TPRT2,  last  time  for  printed  output  (also  end  of 
solution) 

Columns  21-30  DTPRT,  time  interval  for  printed  output 
Columns  31-40  TPLT1.  first  time  for  plotted  output 
Columns  41-50  TPLT2,  last  time  for  plotted  output 
Columns  51-60  DTPLT.  time  interval  for  plotted  output 

C.  Post-processor  DFSTRW 

This  program  creates  deformed  structure  plots  at  specified  times.  Plots  of  the 
entire  warhead  or  “zoomed  in”  plots  of  a  specified  portion  can  be  made.  The  original 
start  tape  for  the  WHAP  run  and  the  displacement  output  tape  for  the  run  are  read 
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by  DFSTRW  from  logical  units  8  and  20.  Additional  input  data  required  is  given 
below. 

1.  ID  card  (6A6) 

Columns  2-5  User’s  code 

Columns  7-36  User’s  name,  phone  and/or  other  identifying  information 

2.  Control  card  (A6.4X.F  10.0,315) 

Columns  1-10  ORIENT,  page  orientation  indicator,  if  ORIENT  = 

“VERTICAL — ”  (or  blank)  plots  will  be  drawn  with 
the  Z-axis  in  the  direction  of  the  larger  page  dimension, 
if  ORIENT  =  “HORIZONTAL”  the  axis  will  be  in  the 
direction  of  shorter  page  dimension 
Columns  11-20  SCALE,  multiplication  factor  for  displacements,  if 
SCALE  =  0.  then  SCALE  =  1. 

Columns  21-25  NFLl,  number  of  first  element  in  block  to  be  plotted, 
if  NELI  =  0  then  NELI  =  1 

Columns  26-30  NELF,  number  of  last  element  in  block  to  be  plotted, 
if  NELF  -  0  then  NELF  =  NEL 

Columns  31-35  ILABEL,  axis  labeling  indicator,  if  ILABEL  >  0  then 
axes  are  drawn  and  labeled 

3.  Plot  area  card  (6F10.0) 

Columns  1-10  RMIN,  R  value  a'  left  end  of  R-axis 

Columns  1 1  -20  RMAX,  R  value  at  right  end  of  R-axis,  if  RMAX  - 

RMIN  =  0.  then  program  will  determine  suitable  values 
Columns  21—30  ZMIN,  Z  value  at  lower  end  of  Z-axis 
Columns  31-40  ZMAX,  Z  value  at  upper  end  of  Z-axis,  if  ZMAX  - 

ZMIN  =  0.  then  program  will  determine  suitable  values 
Columns  41-50  RSTEP,  R  axis  labeling  increment,  if  RSTEP  =  0.  then 
RSTEP  =  RMAX  -  RMIN 

Columns  51-60  ZSTEP,  Z  axis  labeling  increment,  if  ZSTEP  =  0.  then 
ZSTEP  =  ZMAX  -  ZMIN 

4.  Output  segments  card  (1015) 

List  of  output  segments  to  be  plotted,  at  present  only  the  first  and  last 
planes  (planes  1  and  MPLN)  are  plotted  and  must  be  specified. 

5.  Time  cards  (FI 0.0,  one  for  each  plot  desired) 

Columns  1  —  10  T,  time,  if  T  <  0.  then  plot  of  undeformed  structure 
will  be  produced 

D.  Post-processor  CNTRW 

This  program  creates  stress  contour  plots  at  specified  times.  Plots  of  the  entire 
warhead  or  “zoomed  in”  plots  of  a  specified  portion  can  be  made.  The  original  tape 
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tor  the  WHAP  run  and  the  stress  output  tape  for  the  run  are  read  by  C'NTRW  from 
logical  units  8  and  25.  Additional  input  data  required  are  given  below. 

1.  ID  card  (6A6) 

Columns  2-5  User’s  code 

Columns  7-36  User’s  name,  phone  and/or  other  identifying  information 

2.  Control  card  (A6,4X.4I5) 

Columns  1-10  ORIENT,  page  orientation  indicator,  if  ORIENT  = 

“VERTICAL — ”  (or  blank)  plots  will  be  drawn  with 
the  Z-axis  in  the  direction  of  the  larger  page  dimension, 
if  ORIENT  =  “HORIZONTAL”  the  axis  will  be  in  the 
direction  of  shorter  page  dimension 

Columns  11-15  NELI,  number  of  first  element  in  block  to  be  plotted, 
if  NELI  =  0  then  NELI  =  1 

Columns  16—20  NELF,  number  of  last  element  in  block  to  be  plotted, 
if  NELF  -  0  then  NELF  =  NEL 

Columns  21-25  ILABEL,  axis  labeling  indicator,  if  ILABEL  >  0  then 
axes  are  drawn  and  labeled 

Columns  26-30  IBND,  if  IBND  >  0  then  interior  boundaries  between 
different  materials  are  drawn 

3.  Plot  area  card  (6F10.0) 

Columns  1-10  RMIN,  R  value  at  left  end  of  R-axis 

Columns  1 1  -20  RMAX,  R  value  at  right  end  of  R-axis,  if  RMAX  - 

RMIN  =  0.  then  program  will  determine  suitable  values 

Columns  21-30  ZMIN,  Z  value  at  lower  end  of  Z-axis 

Columns  3 1  -40  ZMAX,  Z  value  at  upper  end  of  Z-axis,  if  ZMAX  - 

ZMIN  =  0.  then  program  will  determine  suitable  values 

Columns  41  —  50  RSTEP.  R  axis  labeling  increment,  if  RSTEP  =  0.  then 
RSTEP  =  RMAX  -  RMIN 

Columns  5  1  60  ZSTEP,  Z  axis  labeling  increment,  if  ZSTEP  =  0.  then 
ZSTEP  =  ZMAX  ZMIN 

4.  Output  segments  card  (1015) 

List  of  output  segments  to  be  plotted,  at  present  only  the  first  and  last 

segments  (segments  1  and  MSEC)  are  plotted  and  must  be  specified. 

5.  Contour  description  card  (315, 3F1 0.0.515) 

Columns  1  5  NSICi,  number  indicating  stress  to  be  plotted 
NSIG  Stress  plotted 

1  S  XX 

2  S  YY 

3  S  ZZ 

4  S  XY 

5  S  XZ 
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Si 


6  S  YZ 

7  S  HYDROSTATIC 

8  S-OCTAHEDRAL 


Columns  6-10 

Columns  11-15 

Columns  16-25 
Columns  26-35 

Columns  36  45 
Columns  46—70 


NCNTR.  number  of  contour  levels,  if  NCNTR  =  0  or 
NCNTR  >  10  then  NCNTR  =  10 
J  LABEL,  if  J  LABEL  0  then  contour  numerals  will 
be  connected  with  dashed  lines 

SMIN,  stress  value  corresponding  to  lowest  contour  level 
SMAX,  stress  value  corresponding  to  highest  contour 
level,  if  SMAX  SMIN  -  0.  then  SMIN  and  SMAX  will 
be  determined  by  program  to  span  values  present 
SSTEP,  stress  increment  for  establishing  intermediate 
contour  levels,  overrides  NCNTR 
MATL,  identification  numbers  of  materials  in  which 
contours  are  to  be  drawn  (5  maximum),  if  MATL(l)  = 

0  then  contours  are  drawn  in  all  materials 


6.  Time  cards  (FI 0.0,  one  for  each  plot  desired) 
Columns  1-10  T,  time 
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